home *** CD-ROM | disk | FTP | other *** search
/ Night Owl 6 / Night Owl's Shareware - PDSI-006 - Night Owl Corp (1990).iso / 039a / jpsrc2.zip / JREVDCT.C < prev    next >
C/C++ Source or Header  |  1991-12-02  |  7KB  |  230 lines
  1. /*
  2.  * jrevdct.c
  3.  *
  4.  * Copyright (C) 1991, Thomas G. Lane.
  5.  * This file is part of the Independent JPEG Group's software.
  6.  * For conditions of distribution and use, see the accompanying README file.
  7.  *
  8.  * This file contains the basic inverse-DCT transformation subroutine.
  9.  *
  10.  * This implementation is based on Appendix A.2 of the book
  11.  * "Discrete Cosine Transform---Algorithms, Advantages, Applications"
  12.  * by K.R. Rao and P. Yip  (Academic Press, Inc, London, 1990).
  13.  * It uses scaled fixed-point arithmetic instead of floating point.
  14.  */
  15.  
  16. #include "jinclude.h"
  17.  
  18.  
  19. /* We assume that right shift corresponds to signed division by 2 with
  20.  * rounding towards minus infinity.  This is correct for typical "arithmetic
  21.  * shift" instructions that shift in copies of the sign bit.  But some
  22.  * C compilers implement >> with an unsigned shift.  For these machines you
  23.  * must define RIGHT_SHIFT_IS_UNSIGNED.
  24.  * RIGHT_SHIFT provides a signed right shift of an INT32 quantity.
  25.  * It is only applied with constant shift counts.
  26.  */
  27.  
  28. #ifdef RIGHT_SHIFT_IS_UNSIGNED
  29. #define SHIFT_TEMPS    INT32 shift_temp;
  30. #define RIGHT_SHIFT(x,shft)  \
  31.     ((shift_temp = (x)) < 0 ? \
  32.      (shift_temp >> (shft)) | ((~0) << (32-(shft))) : \
  33.      (shift_temp >> (shft)))
  34. #else
  35. #define SHIFT_TEMPS
  36. #define RIGHT_SHIFT(x,shft)    ((x) >> (shft))
  37. #endif
  38.  
  39.  
  40. /* The poop on this scaling stuff is as follows:
  41.  *
  42.  * We have to do addition and subtraction of the integer inputs, which
  43.  * is no problem, and multiplication by fractional constants, which is
  44.  * a problem to do in integer arithmetic.  We multiply all the constants
  45.  * by DCT_SCALE and convert them to integer constants (thus retaining
  46.  * LG2_DCT_SCALE bits of precision in the constants).  After doing a
  47.  * multiplication we have to divide the product by DCT_SCALE, with proper
  48.  * rounding, to produce the correct output.  The division can be implemented
  49.  * cheaply as a right shift of LG2_DCT_SCALE bits.  The DCT equations also
  50.  * specify an additional division by 2 on the final outputs; this can be
  51.  * folded into the right-shift by shifting one more bit (see UNFIXH).
  52.  *
  53.  * If you are planning to recode this in assembler, you might want to set
  54.  * LG2_DCT_SCALE to 15.  This loses a bit of precision, but then all the
  55.  * multiplications are between 16-bit quantities (given 8-bit JSAMPLEs!)
  56.  * so you could use a signed 16x16=>32 bit multiply instruction instead of
  57.  * full 32x32 multiply.  Unfortunately there's no way to describe such a
  58.  * multiply portably in C, so we've gone for the extra bit of accuracy here.
  59.  */
  60.  
  61. #ifdef EIGHT_BIT_SAMPLES
  62. #define LG2_DCT_SCALE 16
  63. #else
  64. #define LG2_DCT_SCALE 15    /* lose a little precision to avoid overflow */
  65. #endif
  66.  
  67. #define ONE    ((INT32) 1)
  68.  
  69. #define DCT_SCALE (ONE << LG2_DCT_SCALE)
  70.  
  71. /* In some places we shift the inputs left by a couple more bits, */
  72. /* so that they can be added to fractional results without too much */
  73. /* loss of precision. */
  74. #define LG2_OVERSCALE 2
  75. #define OVERSCALE  (ONE << LG2_OVERSCALE)
  76. #define OVERSHIFT(x)  ((x) <<= LG2_OVERSCALE)
  77.  
  78. /* Scale a fractional constant by DCT_SCALE */
  79. #define FIX(x)    ((INT32) ((x) * DCT_SCALE + 0.5))
  80.  
  81. /* Scale a fractional constant by DCT_SCALE/OVERSCALE */
  82. /* Such a constant can be multiplied with an overscaled input */
  83. /* to produce something that's scaled by DCT_SCALE */
  84. #define FIXO(x)  ((INT32) ((x) * DCT_SCALE / OVERSCALE + 0.5))
  85.  
  86. /* Descale and correctly round a value that's scaled by DCT_SCALE */
  87. #define UNFIX(x)   RIGHT_SHIFT((x) + (ONE << (LG2_DCT_SCALE-1)), LG2_DCT_SCALE)
  88.  
  89. /* Same with an additional division by 2, ie, correctly rounded UNFIX(x/2) */
  90. #define UNFIXH(x)  RIGHT_SHIFT((x) + (ONE << LG2_DCT_SCALE), LG2_DCT_SCALE+1)
  91.  
  92. /* Take a value scaled by DCT_SCALE and round to integer scaled by OVERSCALE */
  93. #define UNFIXO(x)  RIGHT_SHIFT((x) + (ONE << (LG2_DCT_SCALE-1-LG2_OVERSCALE)),\
  94.                    LG2_DCT_SCALE-LG2_OVERSCALE)
  95.  
  96. /* Here are the constants we need */
  97. /* SIN_i_j is sine of i*pi/j, scaled by DCT_SCALE */
  98. /* COS_i_j is cosine of i*pi/j, scaled by DCT_SCALE */
  99.  
  100. #define SIN_1_4 FIX(0.707106781)
  101. #define COS_1_4 SIN_1_4
  102.  
  103. #define SIN_1_8 FIX(0.382683432)
  104. #define COS_1_8 FIX(0.923879533)
  105. #define SIN_3_8 COS_1_8
  106. #define COS_3_8 SIN_1_8
  107.  
  108. #define SIN_1_16 FIX(0.195090322)
  109. #define COS_1_16 FIX(0.980785280)
  110. #define SIN_7_16 COS_1_16
  111. #define COS_7_16 SIN_1_16
  112.  
  113. #define SIN_3_16 FIX(0.555570233)
  114. #define COS_3_16 FIX(0.831469612)
  115. #define SIN_5_16 COS_3_16
  116. #define COS_5_16 SIN_3_16
  117.  
  118. /* OSIN_i_j is sine of i*pi/j, scaled by DCT_SCALE/OVERSCALE */
  119. /* OCOS_i_j is cosine of i*pi/j, scaled by DCT_SCALE/OVERSCALE */
  120.  
  121. #define OSIN_1_4 FIXO(0.707106781)
  122. #define OCOS_1_4 OSIN_1_4
  123.  
  124. #define OSIN_1_8 FIXO(0.382683432)
  125. #define OCOS_1_8 FIXO(0.923879533)
  126. #define OSIN_3_8 OCOS_1_8
  127. #define OCOS_3_8 OSIN_1_8
  128.  
  129. #define OSIN_1_16 FIXO(0.195090322)
  130. #define OCOS_1_16 FIXO(0.980785280)
  131. #define OSIN_7_16 OCOS_1_16
  132. #define OCOS_7_16 OSIN_1_16
  133.  
  134. #define OSIN_3_16 FIXO(0.555570233)
  135. #define OCOS_3_16 FIXO(0.831469612)
  136. #define OSIN_5_16 OCOS_3_16
  137. #define OCOS_5_16 OSIN_3_16
  138.  
  139.  
  140. /*
  141.  * Perform a 1-dimensional inverse DCT.
  142.  * Note that this code is specialized to the case DCTSIZE = 8.
  143.  */
  144.  
  145. INLINE
  146. LOCAL void
  147. fast_idct_8 (DCTELEM *in, int stride)
  148. {
  149.   /* many tmps have nonoverlapping lifetime -- flashy register colourers
  150.    * should be able to do this lot very well
  151.    */
  152.   INT32 in0, in1, in2, in3, in4, in5, in6, in7;
  153.   INT32 tmp10, tmp11, tmp12, tmp13;
  154.   INT32 tmp20, tmp21, tmp22, tmp23;
  155.   INT32 tmp30, tmp31;
  156.   INT32 tmp40, tmp41, tmp42, tmp43;
  157.   INT32 tmp50, tmp51, tmp52, tmp53;
  158.   SHIFT_TEMPS
  159.  
  160.   in0 = in[       0];
  161.   in1 = in[stride  ];
  162.   in2 = in[stride*2];
  163.   in3 = in[stride*3];
  164.   in4 = in[stride*4];
  165.   in5 = in[stride*5];
  166.   in6 = in[stride*6];
  167.   in7 = in[stride*7];
  168.  
  169.   /* These values are scaled by DCT_SCALE */
  170.  
  171.   tmp10 = (in0 + in4) * COS_1_4;
  172.   tmp11 = (in0 - in4) * COS_1_4;
  173.   tmp12 = in2 * SIN_1_8 - in6 * COS_1_8;
  174.   tmp13 = in6 * SIN_1_8 + in2 * COS_1_8;
  175.   
  176.   tmp20 = tmp10 + tmp13;
  177.   tmp21 = tmp11 + tmp12;
  178.   tmp22 = tmp11 - tmp12;
  179.   tmp23 = tmp10 - tmp13;
  180.  
  181.   /* These values are scaled by OVERSCALE */
  182.  
  183.   tmp30 = UNFIXO((in3 + in5) * COS_1_4);
  184.   tmp31 = UNFIXO((in3 - in5) * COS_1_4);
  185.  
  186.   OVERSHIFT(in1);
  187.   OVERSHIFT(in7);
  188.  
  189.   tmp40 = in1 + tmp30;
  190.   tmp41 = in7 + tmp31;
  191.   tmp42 = in1 - tmp30;
  192.   tmp43 = in7 - tmp31;
  193.  
  194.   /* And these are scaled by DCT_SCALE */
  195.  
  196.   tmp50 = tmp40 * OCOS_1_16 + tmp41 * OSIN_1_16;
  197.   tmp51 = tmp40 * OSIN_1_16 - tmp41 * OCOS_1_16;
  198.   tmp52 = tmp42 * OCOS_5_16 + tmp43 * OSIN_5_16;
  199.   tmp53 = tmp42 * OSIN_5_16 - tmp43 * OCOS_5_16;
  200.   
  201.   in[       0] = (DCTELEM) UNFIXH(tmp20 + tmp50);
  202.   in[stride  ] = (DCTELEM) UNFIXH(tmp21 + tmp53);
  203.   in[stride*2] = (DCTELEM) UNFIXH(tmp22 + tmp52);
  204.   in[stride*3] = (DCTELEM) UNFIXH(tmp23 + tmp51);
  205.   in[stride*4] = (DCTELEM) UNFIXH(tmp23 - tmp51);
  206.   in[stride*5] = (DCTELEM) UNFIXH(tmp22 - tmp52);
  207.   in[stride*6] = (DCTELEM) UNFIXH(tmp21 - tmp53);
  208.   in[stride*7] = (DCTELEM) UNFIXH(tmp20 - tmp50);
  209. }
  210.  
  211.  
  212. /*
  213.  * Perform the inverse DCT on one block of coefficients.
  214.  *
  215.  * A 2-D IDCT can be done by 1-D IDCT on each row
  216.  * followed by 1-D IDCT on each column.
  217.  */
  218.  
  219. GLOBAL void
  220. j_rev_dct (DCTBLOCK data)
  221. {
  222.   int i;
  223.   
  224.   for (i = 0; i < DCTSIZE; i++)
  225.     fast_idct_8(data+i*DCTSIZE, 1);
  226.   
  227.   for (i = 0; i < DCTSIZE; i++)
  228.     fast_idct_8(data+i, DCTSIZE);
  229. }
  230.